lambda_seq <- seq(1,.305,-.005)^15 # Sequence of candidate penalties to search over
if (rhs.group == 'full' & (country == 'indo')) {
  lambda_seq <- c(seq(14,2,-1), 
                  seq(1,.86,-.01)^15,
                  seq(.85,.635,-.005)^15) # Changed just for an Indo High Run 
  # that was taking forever
  
}

# Initialize index (i will index years)
i = 0
lasso.results <- list()
for (year in start.year:end.year) {
  i = i + 1
  print(year)
  # Fit LASSO
  lasso.results[[i]] <- list()
  dta.past <- dta[dta[,t.var]< year,]
  N <- length(dta.past[,dv])
  cv_results <- foreach (cv=1:cv.runs, 
                         .combine='rbind',
                         .packages = c("glmnet","pROC")) %dopar% {   # Loop over CV runs
     set.seed(52*cv) 
     shuffle <- sample(1:N,
                       N,
                       replace=FALSE)
     pred.prob <- matrix(NA,nrow=N,ncol=length(lambda_seq))
     sampleSplit = c(0,round((N/5*c(1:5))))
     for (f in 1:5) {
       # Fit lasso
       cv.fit <- glmnet(x = data.matrix(dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                    rhs]),
                        y = dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                    dv],
                        alpha = lasso_params$alpha,
                        family = lasso_params$family,
                        lambda = lambda_seq)
       pred.prob[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                 1:length(cv.fit$lambda)] <- predict(cv.fit,
                                                     newx = data.matrix(dta.past[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                                                       rhs]),
                                                                              type = "response")
    }
     cv_results <- c("CV Run"=0,
                     "lambda"=0,
                     "MSE"=10^4)
     # Get error rates for each lambda used
     for (l in seq(10, 10*(round(length(lambda_seq)/10)-1),
                   10)) {
       if (sum(is.na(pred.prob[,l]))==0) {
         mse <- mean((dta.past[,dv]-pred.prob[,l])^2)
         if (mse<cv_results["MSE"]) {
           cv_results <- c("CV Run"=cv,
                           "lambda"=lambda_seq[l],
                           "MSE" = mse)
           best.l <- l
         }
       }
     }
     for (l in seq(best.l-9,best.l+10,1)) {
       if (sum(is.na(pred.prob[,l]))==0) {
         mse <- mean((dta.past[,dv]-pred.prob[,l])^2)
         if (mse<cv_results["MSE"]) {
           cv_results <- c("CV Run"=cv,
                           "lambda"=lambda_seq[l],
                           "MSE" = mse)
           best.l <- l
         }
       }
     }
     cv_results<- c(cv_results,
                    "DT" = as.numeric(mostaccDT(y=dta.past[,dv],
                                                yhat=pred.prob[,best.l],
                                                J=1000)),
                    "DTsens" = as.numeric(bestDT(y=dta.past[,dv],
                                                yhat=pred.prob[,best.l],
                                                J=1000))
                    )
     return(cv_results)
  } 
  hours <- floor((proc.time()[3]-start.time)/3600)
  mins <- floor((proc.time()[3]-start.time)/60) - hours*60
  secs <- floor((proc.time()[3]-start.time)) - hours*3600 - mins*60
  
  print(paste(hours,
              "h",
              mins,
              "m",
              secs,
              "s elapsed",
              sep = ""))
  if (cv.runs > 1) {
    optim.lambda <- mean(cv_results[,"lambda"])
  } else {
    optim.lambda <- cv_results["lambda"]
  }
  
  lasso.results[[i]]$mod <- glmnet(x = data.matrix(dta[dta[,t.var]< year,
                                                          rhs]),
                                  y = dta[dta[,t.var]< year,
                                                         dv],
                                  alpha = lasso_params$alpha,
                                  family = lasso_params$family,
                                  lambda = seq(10.9,1,-.1)^6*optim.lambda
  )
  opt.lam <- 100
  while (is.na(lasso.results[[i]]$mod$lambda[opt.lam])) {
    opt.lam = opt.lam - 1
  }
  lasso.results[[i]]$lambda <- optim.lambda
  if (cv.runs > 1 ) {
    lasso.results[[i]]$dt <- mean(cv_results[,"DT"]) 
    lasso.results[[i]]$dtsens <- mean(cv_results[,"DTsens"])
  } else {
    lasso.results[[i]]$dt <- cv_results["DT"] 
    lasso.results[[i]]$dtsens <- cv_results["DTsens"]
  }
 
  lasso.results[[i]]$fit.oos <- predict(lasso.results[[i]]$mod,
                                       newx = data.matrix(dta[dta[,t.var] == year,
                                                                 rhs]),
                                       type = "response")[,opt.lam]
  lasso.results[[i]]$actual.oos <- as.matrix(dta[dta[,t.var] == year,
                                                      dv])
  lasso.results[[i]]$param_dist <- cv_results
  print(year)
}
save(lasso.results,file=paste(modeldir,"/",
                              country,
                              "_lasso_",
                              v,
                              "_",
                              rhs.group,
                              "_mse",
                              ".RData",
                              sep = ""))